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FOREWORD 
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ter, Hampton, Virginia. The work reported was performed under contract NAS 1-19000. Dr. 
Allan R. Wieting was the NASA Technical Monitor for this task. 

SUMMARY 

This paper presents a nonequilibrium flow solver, implementation of the algorithm on 
unstructured meshes, and application to hypersonic flow past a blunt bodies. Air is modeled as a 
mixture of five chemical species, namely 02, N2, 0, NO, and N, and having two temperatures, 
namely translational and vibrational. The solution algorithm is a cell-centered, point-implicit 
upwind scheme that employs Roe’s flux difference splitting technique. Implementation of this 
algorithm on unstructured meshes is described. The computer code is applied to solve Mach 15 
flow with and without a Type IV shock interference on a cylindrical body of 2.54 mm (0. 1 inch) 
radius representing a cowl lip. Adaptively generated unstructured meshes are employed, and the 
meshes are refined several times until the solution exhibits detailed flow features and surface pres- 
sure and heat flux distributions. Effects of a catalytic wall on surface heat flux distribution are 
studied. For the Mach 15 Type IV shock interference flow, present results showed a peak heat 
flux of 544 MW/m 2 for a fully catalytic wall and 43 1 MW/m 2 for a noncatalytic wall. Some of 
the results are compared with available computational data. 

INTRODUCTION 

The study of hypersonic flows has gained momentum with the advent of concepts like the 
National Aerospace Plane (NASP) and similar transatmospheric vehicles. Under the very high 
velocity and temperature conditions experienced by hypersonic vehicles, departure from chemical 
and thermal equilibrium occurs. Properties of air change dramatically as new chemical species 
are produced at the expense of others. The simple one temperature model used to describe the 
energy of air becomes inapplicable, and it becomes necessary to consider one or more additional 
temperatures (corresponding to vibrational and electronic energies). Determination of aerother- 
mal loads on blunt bodies in such an environment is of great importance, and forms the subject of 
the present study. 



In high speed flows, any adjustment of chemical composition or thermodynamic equilib- 
rium to a change in local environment requires certain time. This is because the redistribution of 
chemical species and internal energies require certain number of molecular collisions, and hence a 
certain characteristic time. Chemical nonequilibrium occurs when the characteristic time for the 
chemical reactions to reach local equilibrium is of the same order as the as the characteristic time 
of the fluid flow. Similarly, thermal nonequilibrium occurs when the characteristic time for trans- 
lational and various internal energy modes to reach local equilibrium is of the same order as the as 
the characteristic time of the fluid flow. Since the chemical and thermal changes are the results of 
collisions between the constituent particles, nonequilibrium effects prevail in high-speed flows in 
low-density air. 

In chemical nonequilibrium flows the mass conservation equation is applied to each of the 
constituent species in the gas mixture. Therefore, the overall mass conservation equation is 
replaced by as many species conservation equations as the number of chemical species consid- 
ered. The assumption of thermal nonequilibrium introduces additional energy conservation equa- 
tions - one for every additional energy mode. Thus, the number of governing equations for 
nonequilibrium flow is much bigger compared to those for perfect gas flow. A complete set of 
governing equations for nonequilibrium flow may be found in Gnoffo, et al. (Ref. 1) and Lee 
(Ref. 2). 

Analysis of nonequilibrium flow is rather complex because (1) the number of equations to 
be solved is much larger than the Navier-Stokes equations, and (2) there are additional terms like 
the species production, mass diffusion, and vibrational energy relaxation, etc., that appear in the 
governing equations. In a typical flight of the NASP flying at Mach 15, ionization is not expected 
to occur, and a 5-species air is adequate for the analysis (see Gupta, et al, Fig. 1, Ref. 3). Since the 
rotational characteristic temperatures for the constituent species (namely 02, N2, O, NO, and N) 
are small, the translational and rotational energy modes are assumed to be in equilibrium, whereas 
the vibrational energy mode is assumed to be in nonequilibrium. We have simplified the thermo- 
dynamic model by assuming a harmonic oscillator to describe the vibrational energy. Ionic spe- 
cies and electrons are not considered. This simplifies the set of governing equations by elimi- 
nating the equation governing electron and electronic excitation energy. In the present analysis, 
we have taken the complete set of governing equations from Gnoffo, et al. (Ref. 1), and simplified 
them for a five-species two-temperature air model. 

A cell-centered, point-implicit, upwind scheme using Roe’s flux difference splitting tech- 
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nique is employed for the inviscid terms of the governing equations. An existing set of computer 
codes, Langley Adaptive Remeshing Code and Havifir Stokes Solver (LARCNESS, see Ref. 4), 
is modified for nonequilibrium flow of air with five chemical species, namely 02, N2, 0, NO, and 
N, and two temperatures, namely translational-rotational and vibrational temperatures. The 
solver is first-order accurate in space and time. Adaptively generated unstructured meshes con- 
sisting of triangular and quadrilateral elements are employed to discretize the computational 
domain; the adaptive remeshing strategy places small elements in regions of large second deriva- 
tives of selected flow variables, and facilitates effective capturing of sharp flow features. 

When an oblique shock wave impinges on a bow shock in front of a blunt body like the 
leading edge of an engine cowl, a complex flow pattern results. Edney (Ref. 5) studied these flow 
patterns extensively, and classified them into six types. Of these, the Type IV shock interference 
(see Fig. 1), is characterized by two triple points, two shear layers, and a supersonic jet (in 
between the shear layers) that undergoes repeated compressions and expansions before impinging 
on the body. Impingement of this supersonic jet causes intense heating at the stagnation point. 
The knowledge of thermal loads forms an important input in the design of such leading edges. 

In the present paper a Type IV shock interference in Mach 15 flow past a blunt body of 
2.54 mm radius is solved. This body is typical of the engine cowl lip of the NASR The effects of 
a catalytic wall on the heat flux distribution on the wall are studied. The results are compared 
with other numerical results. 

Symbols: 

aij constants that appear in computing thermal conductivity 

a t j constants that appear in the curve fits for c p t i 

A Avogadro constant, 6.022 045 E+23 per mol 

c, average speed of molecules of species i, m/s 

c speed of sound, m/s 

Cf = f/p oa uJ‘, skin friction coefficient 

Cp t i specific heat at constant pressure per mole of species i for translational energy, 

J/mol-K 
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c v r i specific heat at constant volume per mole of species i for rotational energy, 

J/mol-K 

c v l j specific heat at constant volume per mole of species i for translational energy, 

J/mol-K 

c v v i specific heat at constant volume per mole of species i for vibrational energy, 

J/mol-K 

D, effective diffusion coefficient for species i 

Di average vibrational energy of molecule i, that is created or destroyed, J/kg 

d distance from the body surface in radial direction, (negative outwards), m 

E =e + (i? + v 2 )/!, total energy (sum of internal and kinetic energies) per unit mass 

of gas mixture, J/kg 

e =e t + e r + internal energy per unit mass of the mixture, J/kg 

e r rotational energy per unit mass of the gas mixture, J/kg 

e t translational energy per unit mass of the gas mixture, J/kg 

e v vibrational energy per unit mass of the gas mixture, J/kg 

e r i rotational energy per mole of species i, J/mol 

e t i translational energy per mole of species i, J/mol 

e v i vibrational energy per mole of species i, J/mol 

e* v i vibrational energy at translational temperature per mole of species /, J/mol 

/ skin friction, N 

G, Gibbs energy for species i, J/mol-K 

H =E + P/p, enthalpy of the gas mixture, J/kg 

h , enthalpy per mole of species i, J/mol 

k Boltzmann constant, 1.380 622 E-23 J/K 

kf,v kb.i forward and backward reaction rate constants for reaction j 

M molecular weight of the mixture, kg/mol 
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T 
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U, V 

Moo 

x i 

y, 

a, 

P 

4 ?- 

Mfi 

T], 

Vr 

Tlv 


molecular weight of species i, kg/mol 
reduced molecular weights of species i and j 
number density of species i, per mol 
pressure, N/m 2 

pitot pressure in perfect gas, N/m 

surface heat flux, MW/m 2 , also 

=('m 2 + v 2 )/?, dynamic pressure, N/m 2 

universal gas constant, 8.314 J/mol-K 

entropy of the gas mixture, J/kg-K 

entropy of species i, J/mol-K 

translational temperature, K 

vibrational temperature, K 

wall temperature, K 

freestream temperature, K 

velocity components in x- and y-directions, m/s 

freestream velocity, m/s 

mass fraction of species i 

rate of production of species i, kg/m 3 -s. 

mole fraction of species i 

dP 

dPt 

dP 

3(p£) 

modified collision integrals for species i and j, m-s 
heat of formation per mole of species i, J/mol 
thermal conductivity for translational energy, J/m-s-K 
thermal conductivity for rotational energy, J/m-s-K 
thermal conductivity for vibrational energy, J/m-s-K 
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density, kg/m 3 

partial density of species i, kg/m 3 

viscosity of the gas mixture, N-s/m 2 

angular location on the body, deg. (see Fig. 2) 

characteristic vibrational temperature of species i, K 

number of moles of the mixture per unit mass of the gas mixture, mol/kg 

number of moles of species i per unit mass of the mixture, mol/kg 

viscous stresses, N/m 2 

translational- vibrational relaxation times for molecular species, s 
effective collision cross-section for vibrational relaxation, m 2 
collision integrals for species i and j, m 2 


Subscripts: 

w 


ij 

r 

t 

v 


wall quantity 
Freestream value 
indices for chemical species 
rotational mode 
translational mode 
vibrational mode 


GOVERNING EQUATIONS 


The governing equations for two-dimensional flow in chemical and thermal nonequilib- 
rium written in conservation form are as follows: 

Species Conservation : 


3p, dP, v 

Bt Bx By 



( 1 ) 
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Momentum Conservation: 


9p« , 3(P + p« 2 ) . dpuv _ 3 x . 3 x 

ir f — ‘ + 5?<V 

3pv 3puv 3(P + pv 2 ) _ 3 . 9 . 

3/ + 3jc + dy - 51 ( V + 


Total Energy Conservation: 


3p E 3p Hu 3p Hv 3 , . 3 , * 

dF + 3 r + 3T " 51 ( “ T ~ + v V + 5y (wXjt > + VT - } 


jy' 


3? 

+ l^l 

Vibrational Energy Conservation: 


'3jp v 3* / 5y'-^ , 3;y dy / 

5 


3y 


V 1 VV^'j 3 [ v hjD.dYi 
2* W t dx 1 5^1 w, dy 

i m 1 7 v i « 1 


(2) 


(3) 


3pe 3pe u dpe v 

- 4 - 4 - _ 

3r 3* 


pZ 

^ i a mo/ 


3y 

K.PpY? 

W, dx 




dy 

KPpYi 

*i *yj 

+ X w (e *‘ -e v.i) /x « + X *A 


Pi 

i » mo/ 


— ,w, 

i * mo/ * 


i ■ mo/ 


(4) 


These equations are simplified form of a complete set of governing equations given in Ref. 1, with 
minor changes in notation. A detailed description of each of the terms in these equations may be 
found in Ref. 1 . 


The reacting gas mixture is assumed to consist of five species, namely 0 2 , N 2 , 0, NO, and 
N; hence there are five components in Eqn. (1) — one for each species. The density, p, of the gas 
mixture is the sum of the partial densities of the five species, i.e., 


5 

p ■ Ip, (5) 

i- 1 

The species production terms, w jt i = l 5 , are in general, non-zero; however, the sum 

£w f = 0 . Therefore, if the diffusion terms are neglected, and all the species conservation equa- 
tions represented by Eqn. (1) are added, then we obtain the familiar overall mass conservation 
equation. 

The gas mixture is assumed to behave as a perfect gas; hence the pressure, P , is given by 
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the equation for a thermally perfect gas, i.e., 

s 

P = pRT^tj, = p RTo (6) 

I « 1 

where o = , and o,, i = 1, .... 5 are the mole numbers of the species i per unit mass of the 

mixture. Note that the molecular weight of the gas mixture, M , is equal to l/o . The value of the 
universal gas constant, R, is 8.3 14 J/mol-K. The above equation can also be written as follows: 

5 

p = (7) 

i-i 

where A# ( , i = l, .... 5 are the molecular weights of species i. 

THERMODYNAMIC MODEL 


The total energy, E , is the sum of the internal energy and the kinetic energy. The internal 
energy is the sum of translational energy, e, , rotational energy, e r , and vibrational energy, e v . 
These energies are given by the following relations: 


( T 


S t = = X 


jc vli (r)dr+Ah fi 


i - 1 


\0 


5 5 T 

= y £o i jc„' l ir)dr 

i* * I i » 1 o 


( 8 ) 


e v ~ - X°-J c v.v.<( r )4r 

i - 1 1-1 o 

The reference temperature is taken as zero Kelvin. The heat of formation, A h f ( for 0 2 and N 2 are 
zero. The values of A h f , for O, NO, and N at the reference temperature taken from Prabhu and 
Erickson (Ref. 6) are 246,783 J/mol, 90,671 J/mol, and 470,816 J/mol, respectively. The heat 
capacities for the three energy modes are given by the following relations: 


c v>u (7) = 2.5R 

Cy.r.iCO = R 


v, , ( T v ) = R 


exp(e v /T v ) 
[exp{Q v i /T v ) - 1] 



(9) 


Note that the translational and rotational energy modes are assumed to be fully excited, and hence 
the heat capacities for these modes are independent of temperature. The vibrational energy mode 
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is assumed to be not fully excited, and hence the vibrational heat capacity is a function of the 
vibrational temperature. The above expression for c v v i is from Vincenti and Kruger (Ref. 7), and 
is the result of the assumption that the molecules may be regarded as harmonic oscillators. The 
characteristic temperature for vibration, 0 V appearing in Eqn. (9) is 2,270 K for O 2 , 3,390 K for 
N 2 , and 2,740 K for NO. Note that when the vibrational energy mode becomes fully excited, i.e., 
when r v » 0 V , c v v i = R. 

The total enthalpy is given by 

H = E + P/p (10) 

In the present computational scheme, the partial derivatives of pressure with respect to the 
conserved variables are required. These derivatives are evaluated as follows. 

First, the equation of state given by Eqn. (6) is rewritten as follow: 

P = TLR i p i (11) 

Note that /?, = R/M i . Equation (11) may be written in differential form as 

dP = RI. (P/A/,.) dT+ TLR i dp i (12) 

Next, an expression for dT can be derived from the definition of total energy, and is given 
by 


dT = \d{pE) -Edp-de v + e v dp-'Le ti {dp i -x i dp)/M i + (u + v 2 ) dp + udpu + v4pv]/ (c v , r p) (13) 
where c v tr = ^a,c v n . . Combining equations (12) and (13), we obtain the following expressions 
for the partial derivatives of pressure: 


dP 

acp E) 



|£« P [(«* 2 + v 2 )/2-e,J+*r/A/ j 


dP 

9(pv) 


-Pv 


= + RT/M i + P (u + v 2 ) /2 

= a ( + p(« 2 + v 2 ) /2 


(14) 


where a, = RT/M-$e lr i . These partial derivatives of P will be used in the solution algorithm. 

Speed of sound is another variable required in the present analysis. An expression for the 
speed of sound can be obtained starting from the first principles. In problems with nonequilib- 
rium flows, the appropriate speed of sound is the frozen speed of sound (see Ref. 7) defined by 
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c 2 = dP 


dp 


(15) 


S,9 


where 5 is the entropy of the gas mixture and $ are the nonequilibrium variables, namely the 
mole numbers of the species and the vibrational energy. For a mixture of gases, the entropy is 
given by 


S = ^sfl,-Ra[ln (p/p 0 ) -^/n(o/o ( )] 
and the entropy, s { , of each of the species is given by 


dh. 


S, = 


(16) 


(17) 


If the vibrational energy is held constant, then we can write 

- c p.J-f (18) 

Differentiating Eqn. (16) with constant a, , and substituting for ds t from Eqn. (18), we obtain 

dS = 'L c p.,r.^Y- R<J T (19) 

Also, upon differentiating the equation of state with constant o i yields 


dP _d£ dT 
P p T 

Now setting dS = 0 , eliminating dT between Eqns. (19) and (20), and using the relation 

= R ~ c v,tr,i> we find 


( 20 ) 


dP _ dp Ra dP _ do P dP 

P ~ p S W*/ = p p +1p 

Hence 

cJ = (DL, = (P+1,,,/p 

Next, consider the expression 


( 21 ) 


( 22 ) 


+ p/»-p* v = ^(RT/M r ^e, r>i )x i + p*-pe v 
= RTa - Pe fr + pA-P* v 
= P/p - Pe, r + P (*„ + e v + P/p ) ~P« V 
= (P+DF/p 

Therefore, the speed of sound in nonequilibrium flow can be expressed as 


( 23 ) 
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( 24 ) 


c = 7X a « :t . + p (/,_e v ) 

TRANSPORT PROPERTIES 

A procedure for computation of the transport properties of a multicomponent gas mixture 
has been described by Gnoffo, et al. (Ref. 1), which in turn is based on Lee’s formulation (Ref. 2). 
These computations require evaluation of two collision integrals. These collision integrals 
depend only on the translational temperature, 7\ for atomic and molecular particles, (electrons are 
absent in the present analysis), and are evaluated as curve fits using the tabulated data also from 
Ref. 1. We follow the procedure of Gnoffo, et al. (Ref. 1) for the computation of viscosity and 
thermal conductivity of the gas mixture, and binary diffusion coefficients of the 5 species. This 
procedure is summarized here. 

Values of the two collision integrals *\ k = l, 2 , at temperatures of 2,000 K and 4,000 
K are available in tabular form in Ref. 1 . Assuming a linear variation for log [tcQ**' **] (T) with 
log(T) , the following expression is obtained: 

logOn,;^] (T) = log [xQ^] (2000) + slope x In (772000) (25) 

where 

slope = [log [nQ^ k) ] (4000) - log [nQ^* ) ] (2000)]/fn (2) (26) 

The modified collision integrals used in the evaluation of the transport properties are defined as 
follows: 


A (1 ) m 8f 2m i,;1 1/2 ^^(l,l) 

A u M = sL^rJ *°w 
.( 2 )._ 16 [ 2 ^ 1 1/2 ( 2 , 2 ) 


where m , is the reduced molecular weights of species i and j, and is defined by 

li y ^ 

m u = M.M/ (M i + M.) 


Viscosity: 

Viscosity of the gas mixture is given by 


(27) 

(28) 


(29) 
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( 30 ) 


X m i o i 

- — 5 — 

'" 1 2a a S >(7 > 

j - 1 

Thermal Conductivity: 

The translational energy mode is assumed fully excited; hence the translational energy 
thermal conductivity is given by 


n. = 

where the constant a, . is defined as 

*•/ 


, M/M,) [ 0.45 - 2.54 (M/M} ] 

. (./ * + ? 

( 1 + A/,/ M/ 2 

The rotational energy mode is also assumed to be fully excited; hence the rotational energy ther- 
mal conductivity is given by 


= ^ (33) 

i * mol TT* (1) 

2- C A; (^ 

j m 1 

The frozen thermal conductivity is the sum of r\, , and i\ r . The vibrational energy mode, however, 
is assumed to be partially excited; hence the vibrational thermal conductivity is computed follow- 
ing the procedure given in Gupta, et al. (Ref. 3) as follows: 


hv = * X 

i * mol 


(c vvi /R)a. 

3 




(34) 


However, when the vibrational energy mode becomes fully excited, c v v i = R, and n v = which 
is the relation given in Ref. 1. 


Mass Diffusion Coefficient: 

The effective mass diffusion coefficient, D, , of species i in the gas mixture is given by 
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( 35 ) 


„ 0 2 M,(1-0M 

U i “ 5 

Iv\ 

i » < 

where D (J is the binary diffusion coefficient for a pair of particles of species i and j, and is related 
to the modified collision integral as follows: 


D u = 


kT 


P*l'j M 


VIBRATIONAL AND RELAXATION PROCESSES 


(36) 


Millikan and White (Ref. 8) give the following semi-empirical correlation for the vibra- 
tional relaxation time 


= j Xv^ [ A >( 701/3 - ooi5m 'j 4 ) - 18 - 42 1 } / X n > ( 37 ) 

7-1 J i-i 

where nti j is the reduced molecular weight defined in Eqn. (29). The values of the constants Aj 
are 129 for O 2 , 220 for N 2 , and 168 for NO. The relaxation times computed by the above relation 
are valid over a temperature range of 300 to 8,000 K. At temperatures beyond 8,000 K, the above 
relation yields relaxation times much smaller than observed in experiments. For temperatures 
above 8,000 K, Park (Ref. 9) suggests the following relation for the vibrational relaxation time: 

if = l/(*V» ( .e ( .) (38) 

where Q, is the effective collision cross section for vibrational relaxation (following Park, Ref. 9, 
this is set equal to 10‘ 16 cm 2 ), n i is the number density of species i , and ? ( is the average velocity 
of molecules i, given by 

f ( . = [8*7V( Jim,.)] 1/2 (39) 

Blending of the two relations Eqns. (37) and (38) gives the following expression for the vibra- 
tional relaxation time: 


x, = < + < (40) 

Park points out that this expression for vibrational relaxation time is applicable over a much wider 
temperature range. 
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CHEMICAL KINETIC MODEL 


The nonequilibrium chemistry effects involving the five chemical species are modeled by 
17 kinetic steps listed in Table 1. The rate of production of species i is given by the following 
expression: 

17 

= X < v w ( R fj- R bj ( 4l ) 

i-i 

where v. y and v’,^ are the stoichiometric coefficients for the i-th species in the y-th reaction for 
reactants and products, respectively, and R fj and R bj are the forward and backward reaction rates, 
respectively, for the y'-th reaction. The reaction rates are given by the following: 


5 

i ■ 1 
5 

r ».j = vn ( p/ M <> Vt/ 

i« i 


(42) 


where k f j and k b j are the forward and backward reaction rate constants. These rate constants are, 
in general, functions of the translational and vibrational temperatures, see Park (Ref. 10), for 
example. In addition, there is a chemical-vibrational coupling. However, in the present analysis, 
the reaction rate constants are assumed to be functions of translational temperature only, and are 
expressed as follows: 


k f J = AjT^exp (,-dj/T) (43) 

The constants A jt n, , and dj for the 17 kinetic steps are obtained from Gnoffo and McCandless 
(Ref. 11), and are also listed in Table 1. The backward reaction rate constants are obtained from 
the relation 

*',i = h/Ki ( 44 ) 

K e j is the equilibrium constant for the kinetic step y. These equilibrium constants for the 17 
kinetic steps are determined using the following relation: 

K tj = exp (-AGj/RT) (45) 

where A G ; is the net change in the Gibbs energy in reaction y, and is expressed as 

AG , = X (v VV G ‘ (46) 
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where G, is the standard state Gibbs energy of formation of the i-th species. This may be com- 
puted (by noting that G i = A, - Ts t ) by the following expression: 


G f = RT 


a i. 1 ( 1 - ,0 8 ( ' T) ) + a i, 6 /T ~ a i. 7 - X a i./ ,1/ - / (/ - *) 


i-2 


(47) 


Note that contributions from only the translational temperature to A, and s t are considered. The 
constants a i jt i = l, ..., 5 J = l, .... 7 are the constants that appear in the curve fits for the heat 
capacity of the five species, and are taken from Prabhu and Erickson (Ref. 6). 


METHOD OF SOLUTION 


The method of solution employed is this paper is based on the method described by 
Thareja, et al. (Ref. 4) which is an unstructured mesh implementation of a cell centered upwind 
scheme using Roe’s flux difference splitting technique for the inviscid terms of the governing 
equations for perfect gas. The procedure has also been extended to equilibrium flows by Prabhu, 
et al. (Ref. 12). We briefly describe the procedure as it is applied to the equations governing non- 
equilibrium flows. 

The governing equations are written in the following compact form: 


9V.9F.BG _ w 


(48) 


Note that the transport terms are not included here; they will be treated separately, and 


U T = {pj.pK, pv, p£, pc v } 

F T = { p ( u, P + p u, pwv, p uH, p ue v ] 

T 2 W 

G = {p ( v, p«v,P + pv ,pvW, pve v } 

W T = {vv,.. 0,0,0, £ p, (<,. ~e vJ )/x i+ £ h-A) 

i * mol I * mol 

Also note that the source like terms on the right hand side of the vibrational energy conservation 
equation have been included in the vector W. These terms are: (1) energy relaxation between 
vibrational and translational modes, and (2) vibrational energy lost/gained due to molecular 
depletion/production within a cell. 

The computational domain is subdivided into triangular and/or quadrilateral elements. 
The unknowns, V , are assumed to be constant over each element, and the governing equations are 
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integrated over each element. 




dy 

The time-derivative is linearized, and the above equation is rewritten as follows: 


(50) 




(51) 


Indices m and m + 1 refer to two consecutive time levels. Upon replacing the area integral of the 
derivatives with the appropriate contour integrals, we obtain 


^a, + jF^* l dr=w m * 1 a t (52) 

The quantity F n is the component of the inviscid flux normal to a side. The unknowns are discon- 
tinuous across each side. The flux across such a face is given by 

= [i r "t 1 +^* , -W(tC +1 - £ C +, J]/2 (53) 

where the subscript e refers to the element under consideration and the subscript r refers to the 
element sharing the common side. The quantities F^ t and F„ r are the fluxes normal to the side 
under consideration. The matrix |A| is defined as 

\A\ = *|A|* _1 (54) 

where R and R~' are the right and the left eigenvector matrices of the flux Jacobian matrix, A . 
The matrix |A| is the diagonal matrix of the eigenvalues of A , with all negative values replaced by 
corresponding positive values. This introduces upwind bias. The matrix A must be evaluated at 
each element interface, and should exhibit Roe’s property U (see Roe, Ref. 13). Upon substitut- 
ing Eqn. (53) in Eqn. (52), replacing the contour integral by a summation over the sides of the ele- 
ment, and rearranging, we find 

X l>« 1 + iC + 1 - K c? + ‘-tC + 1 )] 5 / 2 = (55) 

sides 

where b s is the length of the side. Linearization of the above equation gives the following 
scheme. 


r /+ ^f5>* |8 ’l At/ = + (56) 

L e s J e 3 

The values with an asterisk as the index are those computed using most recent values in the ele- 
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ment and its neighbors. This makes the scheme point implicit. Evaluation of the matrix A at the 
cell interface forms a major part of the computations, and is discussed in the following para- 
graphs. 

THE JACOBIAN AND RELATED MATRICES 


Gnoffo, et al. (Ref. 1) have derived the inviscid flux Jacobian and eigenvector matrices for 
nonequilibrium flows in 3-D. They include ionic species and electrons in their chemistry model. 
We adopt their formulation simplifying it for 2-D flows with only five chemical species. The 
component of the inviscid flux normal to a side of an element is given by 


F n = Fn x + Gn y (57) 

where n x and n y are the x- and y-components of the unit outward normal to the side. The inviscid 
flux Jacobian defined by A = 3 F n /BU can be derived, and is given by 



■ i\ r x t u 

X i n x 

X ‘ n y 

0 

0 


(a y + P?) n x -Uu 

(1-p) un x + U 

(1-P) vn x - V 

K 

-P"x 

[A] = 

(dj + Pq)n y - Uv 

(1-P)un y +V 

( 1 - P) vn y + U 

P«, 

-p« y 


(a^^q-H)U 

Hn x - pt/« 

Hn y - Pt/v 

d + p)t/ 

-pi/ 


. ~ Ue v 

e v n x 

e v n y 

0 

u 


(58) 


The Kronecker delta function, S ij has the usual meaning, i.e., 8 tJ = 1 if i = j, 5 (J = 0 otherwise. In 
Eqn. (58), the index i=l,...,5 refers to the rows and the chemical species, and the index y'=l,...,5 
refers to the columns and the chemical species. Also U = un x + vn y , V = vn x - un y , and 
q = (h 2 + v 2 ) / 2 . The matrix A is a square matrix of order 9. The eigenvalues of A are 


A = [U ( repeated 6 times) , (U + c) , (1/ - c) , U\ (59) 

where c is the speed of sound given by c = + $(H-q-e v ) . Note that only three of these 

values, namely u,U + c, and U-c, are distinct. The right eigenvector matrix of A corresponding 
to the eigenvalues given by Eqn. (59) is 



h, 

0 

V 2 

x /2 | 

IR] - 

2 

U 


(u + cn x )/2 

(u-cn x )/2 | 

V 

“"x 

(v + cn y ) /2 

(v + cn y )/2 ( 

C 

(P <?-<*,.) 

V 

(H + c 10/2 

(H-cU)/ 2 


0 

0 

V 2 

e/2 

The inverse of the matrix R given by R 1 

is also required, and is given by 


17 



[*-'] = 


( 61 ) 


c 2 8y-x,.(a ; + p 9 ) 

-P«*i 

-pvx, 

-P x i 

P x i 

-V 

~ H y 

n x 

0 

0 

0 lj + P q+Uc 

cn x - J5« 

cn y - Pv 

P 

-P 

a . + P q-Uc 

-cn x -$u 

-c«j,-Pv 

P 

-P 

-e v (Cty+P<?) 

-P«<\, 

-Pv« v 

-P«v 

c + Pe, 


It can be verified that RAR' 1 = A, and RR~' = /, where / is the identity matrix. Note that the 
matrices A, R, and J?" 1 must be evaluated at the so called Roe-averaged state. The variables at 
the Roe-averaged state are designated by an overbar, and are determined as follows: 

One of the properties that the matrix A = A (£/) must satisfy is 


A(F„)=AA(U) (62) 

where A ( ) = ( ) R - ( ) t . Equation (62) has nine components corresponding to the nine compo- 
nents of A (F n ) . The first five components of these are satisfied if we define 


where 


xi - a (*,) L + b (x t ) R i = l 5 

u = au L + bu R v = av L + bv R 


(63) 


a = JPl / (Jp'k + b = + JPl) (64) 

Note that a + b = 1 , and that 


X*. = «X(*<) t + *2 (*«)* = « + *= 1 (65) 

This is a necessary condition if the x, ’s are to retain the property of mass fraction. The sixth and 
the seventh components of Eqn. (62) are satisfied provided the following condition is satisfied. 


s 

A(P) = ^a,A(p,.) +pA(pe„) 

i*l 


Finally, the eighth and the ninth components are satisfied if we define 


( 66 ) 


H = aH,+ bH„ 

(67) 

«V = a(e v ) L + b(e v ) R 

It now remains to determine a, and p . The method of dete rmining these is not obvious, 
since the only condition they have to satisfy is Eqn. (66). In all the available literature, there is no 
rigorous method for the determination of these quantities. Methods based on either some approx- 
imations (see Liu and Vinocur, Ref. 14) or on computational experience (see Gnoffo, Ref. 15) 
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have been suggested. We describe here a simple, yet accurate method to determine these quanti- 
ties. The basis for the present method is the requirement that a,- and p together with the other 
Roe-average variable form a consistent set of thermodynamic variable. 

Since p depends only on the mass fractions, x, , we require that p be related to x, by the 
same functional relation. Hence, P can be determined from the known values of x , , i.e.. 


r J 


P = 


Li = l 


r 5 


Xv Af, / 


( 68 ) 


where ^ = 2.5 for molecular species, and 1.5 for atomic species. Similarly, we require that the 
relation between a, and other variables be the same as the relation between a, and other vari- 
ables. Since a, = RT/M r $e tr i , we write 


a, = RT/M r $e tr j (69) 

Note that c, r , ( the is sum of translational, rotational, and formation energies, and is related to the 
correspond temperature T by the following relation: 


e, r i = (2.5RT + Ah f ,)/Af for molecules 

/.<' < (70) 

= ( 1 .5RT + A h f ,) /M t for atoms 

Therefore, a„ i = l 5 are linear functions of f. Upon using these relations in Eqn. (66), and 

solving the resulting equation for f, we obtain the following: 

5 

A (P) - pA (pe fr ) + P X, A^>/ ,A ( p ; ) / Af, 

T= j ^ (71) 

^Rd-ffrAipJ/Mi 

i* 1 

where, as before,^- = 2.5 for molecular species, and 1.5 for atomic species. All the terms on the 
right hand side of the above equation are known; hence, T can be readily computed. Once f is 
determined, e, r< < can be determined using Eqn. (70). Next, Eqn. (69) provides the values of a, . 

Thus, a „ i = l 5 and p are computed such that the necessary condition Eqn. (66) is satisfied 

exactly. 

Before proceeding further, it is necessary to show that the temperature, f, computed by 
the above equation is always positive. We demonstrate this in the following manner: 

We note that A (P) = (P) R -(P) L , A(p,) = (p,),,- (p,) t , and A(pe, r ) = (pe lr ) R - (pe tr ) L . 
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Substituting these relations, and the expression for p from Eqn. (68) in Eqn. (71), and upon con- 
siderable simplification, we find 

T=aT L + bT R (72) 

This is a remarkably simple result, and confirms that T is positive, and bounded by T L and T R . 
The speed of sound can now be determined using the following relation: 

c = + (73) 

where q = [u 2 + v 2 }/2. Since all the quantities on the r. h. s. of Eqn. (73) are known, c can be 
computed. 

The value of c computed by Eqn. (73) must be positive. To verify this, we start with the 
definition of H 

H = aH L + bH R = ah L + bh R + q + A (74) 

where, it can be shown that 

A = ab[{u L -u R ) 2 + (y L -v R ) T ]/2'2. 0 (75) 

Next, we substitute h = e lr + e v + P/p in Eqn. (74), and obtain the following: 

H = a[(e lr ) L + (<? v ) t + ( p/ P) J +M («„•)* + (<V) * + ( f/ P)*l +9 + 4 
= e lr +e v + (P/p) +q + A 

where e lr = a(e tr ) L + b(e tr ) R and (P/p) = a(P/p) L + b(P/p) R . Next, we consider the definition 
of c and write 


c 2 = X a + P (H - q - e v ) (77) 

Upon substituting for a, and H from Eqns. (69) and (76), respectively, and simplifying, we find 

C 2 = P(e+ (P7p) +A) (78) 

where e = at L + bt R , and e = RT^/p/W , . All the quantities on the right hand side of Eqn. (78) 
are positive, hence c is positive. 

Thus, all the unknowns in the Roe-averaged state required to compute the matrices A , R , 
and p ' are determined. No approximations were made determining a, p , and other variables at 

the Roe-averaged state. These variables namely, jc„ i = 1 5, «, v, H, e„ a„ i = 1, .... 5 , p and c 

form a consistent set of thermodynamic variables, and satisfy the necessary condition Eqn. (66) 
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exactly. 

TRANSPORT TERMS 

All the transport terms, namely the viscous, heat conduction, and mass diffusion terms in 
the governing equations are treated implicitly in essentially the same way as in Thareja, et al. 
(Ref. 4). Similarly, the species production terms are also treated implicitly. These details are not 
included here. The reader is directed to Ref. 4 for information pertaining to the treatment of trans- 
port terms. 

RESULTS AND DISCUSSIONS 


Mach 15 flow past a 2.54 mm radius cylindrical body with and without shock interference 
were considered. The following examples were solved: 

(1) Undisturbed flow 

(2) Type IV shock interference flow 

In both cases, the freestream conditions were assumed to be as follows: 


Velocity 
Density 
Temperature 
Mach number 
Reynolds number 


4,678 m/s 
0.00922 kg/m 3 
241 K 
15 

7,500 


Since the Reynolds number is low, the flow is assumed to laminar everywhere. For the shock 
interference case, the impinging oblique shock angle is 8.88 deg. The location of the impinging 
shock is assumed to be as shown in Fig. 2. Other flow properties in that stream are as follows: 


Velocity 
Density 
Temperature 
Mach number 
Flow inclination 

Wall temperature for both the cases is assumed to be 81 1 K. 


4,628 m/s 
0.02863 kg/m 3 
473 K 
10.6 

6 deg. (up) 

These conditions as well as the loca- 
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tions of the impinging shock are very close to the conditions used by Carlson and Wilmoth (Refs. 
16 and 17), who studied these problems using Direct Simulation Monte Carlo (DSMC) method. 
This facilitated a direct comparison of present results with those from Carlson and Wilmoth. Both 
fully catalytic and noncatalytic walls were considered in the present work, and the wall catalytic 
effects are evaluated. A fully catalytic wall was assumed to have the property that all the atoms 
that came into contact with it recombined into molecules. The fully catalytic wall was assumed to 
have no effect on nitric oxide. 

Mesh Requirements: 

The adaptive remeshing strategy of LARCNESS yields satisfactory meshes for flow prob- 
lems dominated by shock waves. It provided adequate meshes for the case of undisturbed flow 
(with no impinging shock) past the body. Since the undisturbed flow is symmetric about the x- 
axis, flow over 1/2 of the body (1/4 cylinder) was computed. Computations were made on a 
sequence of meshes, each finer than the previous until the peak heat flux on the body remained 
practically unchanged on two successive meshes. The final mesh used in this case (not shown) 
has a total of 17,884 elements consisting of triangles and quadrilaterals. 

In a Type IV shock-shock interference flow the primary flow features are two triple points, 
two shear layers that originate at these triple points, and a supersonic jet that undergoes repeated 
expansion and compression. The procedure of adaptive remeshing employed in LARCNESS, 
captures the strong shocks, but does not capture the weaker shocks and the shear layer. With the 
adaptive remeshing it was possible to capture the bow shock and the transmitted shock in front of 
the body. The regions covered by shear layers and the terminating normal shock, however, 
needed local mesh enrichments. 

As in the previous example, computations were made on a sequence of meshes, each finer 
than the previous until the peak heat flux on the body remained practically unchanged on two suc- 
cessive meshes. The contour plots of flow variable from an initial mesh provided information for 
adaptive remeshing as well as for refinement. The shear layers were identified from the («/« M ) 
contour plots, and these regions were uniformly refined. Uniform refinement subdivides each tri- 
angular element into four smaller ones. The final mesh used in the present computations has the 
following properties: 
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Total number of nodes 
Total number of elements 

Number of triangular elements 
Number of quadrilateral elements 
Total number of sides in the triangles mesh 
Number of sides with length < 0.00254 mm 
Number of sides with length < 0.00508 mm 


85,293 

155,960 

142,069 

13,891 

213,441 

11,876 

82,453 


In this mesh, a large fraction of the total number of elements is in the region of the two 
shear layers and the supersonic jet. The length of the upper shear layer was measured from the u- 
velocity contour plots; its length was 1.0 mm (approx.). It may be recalled that the radius of the 
cylindrical body is 2.54 mm. Following the procedure of Glass (Ref. 18), the shear layer thick- 
ness was estimated to be 0.076 mm, using the known flow conditions downstream of the bow 
shock and assuming laminar flow. (It may be recalled that the ffeestream Reynolds number based 
on the body radius is 1500; hence it was assumed that the flow is laminar everywhere.) Similarly, 
the thickness of the lower shear layer was estimated to be 0.0254 mm. (Note that these shear 
layer originate at their respective triple points with zero thickness.) Triangular elements with sides 
of the order of 0.0025 mm were considered adequate to resolve flow in shear layers of this thick- 

ness. 

The mesh has 29 layers of quadrilateral elements next to the body. This part of the mesh 
helps in accurate determination of the velocity and temperature gradients in the boundary layer. 
The thickness of the first layer is 2.54x1 0' 6 mm, and the total thickness of the quadrilateral ele- 
ments layer is 0.245 mm. There are 480 points along the circumference of the body, with an aver- 
age spacing of 0.376 deg. between points, and a minimum spacing of about 0.056 deg. near the jet 
impingement region. This spacing captured the steep temperature gradients in the stagnation 


region. 


Example 1 - Undisturbed Flow: 

The results for Mach 15 flow past a 2.54 mm radius cylindrical body are shown in Figs. 3- 
11, for the fully catalytic as well as the noncatalytic walls. The contours of (u/uj, Mach number, 
and translational and vibrational temperatures for the fully catalytic wall case are shown in Figs. 
3(a) - 3(d). These plots display smooth contours of flow variables. The maximum translational 
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temperature in the flowfield is 9,937 K, and the maximum vibrational temperature is 5,25 1 K. 

Contours of 02, 0, N2, and N mole fraction for the fully catalytic wall case are shown in 
Figs. 4(a) - 4(d). These plots indicate that oxygen is almost completely dissociated behind the 
normal part of the bow shock. However, only a small part of nitrogen is dissociated. Figures 5(a) 
& 5(b) and 6(a) & 6(b) show similar plots for the noncatalytic wall. No noticeable differences 
may be seen between these plots and the corresponding plots for the catalytic wall. The maxi- 
mum translational temperature is 9,918 K and maximum vibrational temperature is 5,256 K 
(which are close to the corresponding values for the fully catalytic wall). The wall catalytic 
effects are not felt away from the wall. 

The variation of mole fractions along the center line are shown in Figs. 7(a) and 7(b) for 
the fully catalytic and the noncatalytic walls, respectively. As noted earlier from the contour 
plots, oxygen is almost completely dissociated between the bow shock and the body. However, 
all of the oxygen atoms recombine as the flow reaches the fully catalytic wall. In the case of the 
noncatalytic wall, although the wall is relatively cold (T w = 811 K), only some of the oxygen 
atoms recombine into oxygen molecules. This is the result of chemical nonequilibrium. Between 
the bow shock and the body, there is some nitrogen dissociation; but all of the nitrogen atoms 
recombine for fully catalytic as well as the noncatalytic walls. For the noncatalytic walls, the 
recombination of nitrogen atoms is due to the low wall temperature. Mole fraction of nitric oxide 
is not noticeably affected by the wall catalytic effects. The effect of the wall catalyticity is con- 
fined to a region very close to the wall. The shock stand-off distance in both cases is about 0.26 r. 

Variations of (u/tO, (T/T^), and ( T/T „,) along the centerline (line of symmetry) are 
shown in Figs. 8 and 9, respectively. The u- velocity drops sharply as the flow passes through the 
shock, after which it decrease steadily and reaches zero at the stagnation point. The translational 
temperature rises steeply through the shock, and reaches its peak value immediately behind the 
shock. In contrast, the vibrational temperature rises slowly through the shock. This trend is to be 
expected. Behind the shock, as the flow goes towards the stagnation point, the translational tem- 
perature drops from its peak value and the vibrational temperature rises. This is the consequence 
of thermal nonequilibrium. As the flow approaches the stagnation point, both the temperatures 
drop sharply, and reach the specified value of 81 1 K at the wall. The wall catalytic effects have 
very little effect on velocity and temperature variations along the center line. Very close to the 
wall, however, there are significant differences in the temperature distributions between the cata- 
lytic and the noncatalytic walls which lead to differences in the surface heat fluxes. 
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Distribution of pressure amplification (f/P^) on the body is shown in Figs. 10. It may be 
observed that the pressure distribution on the body is unaffected by the wall catalytic effects. At 
the stagnation point the value of {P/P a) is 1.028 implying that the pressure at the stagnation point 
for this nonequilibrium flow is about 2.8% higher than the corresponding value for perfect gas. 
The losses in the total pressure through the normal shock in front of the body, seem to be smaller 
for nonequilibrium flow compared to the perfect gas case. The value of P ,2 is easily computed, 
and is equal to 186 kPa. Hence the stagnation point pressure is 191 kPa. 

The heat flux distribution on the body for the fully catalytic and noncatalytic walls is 
shown in Fig. 11. The surface heat flux for the fully catalytic wall is significantly higher than for 
the noncatalytic wall. The peak values are 20.4 MW/m 2 for the catalytic wall and 15.9 MW/m 2 
for the noncatalytic wall. Unlike the pressure distribution, the heat flux distribution curves are not 
smooth, particularly near the stagnation point. This behavior could not be explained. The very 
large aspect ratio of the elements near the stagnation point is suspected to be one of the causes of 
this. 

Carlson and Wilmoth (Refs. 16 and 17) solved this problem using a DSMC approach and 
assuming a noncatalytic wall. The stagnation point pressure and heat flux obtained by them are 
200 kPa and 15 MW/m 2 , respectively, which compare very well with the corresponding values of 
191 kPa and 15.9 MW/m 2 obtained by the present computations. 

Example 2- Type IV Shock Interference Flow: 

The results for Type IV shock interference flow at Mach 15 past a cylindrical body of 2.54 
mm radius are shown in Figs. 12 - 20. The mole fraction contours of 02 and O are shown in Fig. 
12, and the mole fraction contours of N2 and N are shown in Fig. 13. These plots are for the fully 
catalytic wall; corresponding plots for the noncatalytic wall are not presented as they were not 
noticeably different from the plots for the fully catalytic wall except very close to the wall. The 
wall catalytic effect does introduce differences in the gas composition close to the wall. This will 
be discussed later in this section. As in the case of undisturbed flow, most of the oxygen is disso- 
ciated behind the normal shock, whereas only a small amount of nitrogen is dissociated. 

Contour plots of (u/iO and Mach number are shown in Fig. 14, and the contour plots of 
translational and vibrational temperature are shown on Fig. 15. These contour plots are for the 
catalytic wall; there are no noticeable differences between these contours plots and the corre- 
sponding plots for the noncatalytic wall. The two triple points, the two shear layers, and the jet 


25 



terminating normal shock may be seen in these figures. The maximum translational temperature 
is 10,082 K, and is reached as the flow passes through the normal part of the bow shock above the 
impinging oblique shock. The maximum vibrational temperature is 5,594 K. The maxim um 
translational temperature is close to the corresponding value for the undisturbed flow (9,937 K). 
This is to be expected, since the maximum temperature occurs as the flow passes through the nor- 
mal shock. 

The pressure contours in the flow field are shown on Figs. 16a and 16b. The two triple 
points, the transmitted shock, and a compression and an expansion of the flow (before it reaches 
the terminating normal shock) may be seen in Fig. 16a. An enlarged view of the pressure con- 
tours near the stagnation point is shown in Fig. 16b. The normal shock that terminates the super- 
sonic jet may be seen clearly in this figure. This normal shock is essentially parallel to the wall, 
and hence is expected to lead to the most severe heating on the body. The peak pressure on the 
body is at 0 = -16.7 deg. The pressure on the body over a one-degree interval (0 = -16 deg. to -17 
deg.) is essentially constant. It may be noted that this region lies right behind the normal shock. 
The pressure contours for the noncatalytic wall are not shown; there is no noticeable differences 
between the pressure contours for the catalytic and noncatalytic wall cases. 

The pressure distributions on the body for the catalytic and the noncatalytic walls are 
shown in Figs. 17a and 17b. The wall catalytic effects do not seem to affect noticeably, either the 
peak pressure or the pressure distribution on the body. Although the pressure distribution appears 
to have a sharp peak, Fig. 17a, the peak value is essentially constant over a one-degree range - 
between 0 = -16 deg. and -17 deg. (See Fig. 17b). As noted earlier, this corresponds to the region 
covered by the terminating normal shock of the supersonic jet. The value of the pressure peak is 
about 13.1 times the pressure at the stagnation point for undisturbed flow of a perfect gas, and 
occurs at 0 = -16.7 deg. on the body (also see the pressure contours, Fig. 16b). With P t2 equal to 
186 kPa, the maximum pressure is 2,440 kPa. 

The heat flux distributions for the catalytic and the noncatalytic walls are shown in Figs. 
18a and 18b. Similar to the pressure distribution, the heat flux distribution also exhibits a sharp 
peak (See Fig. 18a). A closer examination of this peak, however, reveals some details (Fig. 18b). 
Unlike the pressure distribution, the heat flux distribution for the fully catalytic wall is signifi- 
cantly higher than for the noncatalytic wall. Also, the distributions have two peaks - one at 0 = - 
16 deg. and the other at 0 = -17 deg. Recall that this region (between 0 = -16 deg. and -17 deg.) 
over which the pressure is practically constant, corresponds to the region covered by the normal 
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shock terminating the supersonic jet. The probable cause of this particular shape of the heat flux 
distribution with two peaks is discussed in a following paragraph. The peak value is 544 MW/m 2 
for the catalytic wall and 43 1 MW/m 2 for the noncatalytic wall. At the fully catalytic wall, all the 
atoms recombine into molecules. Atomic recombination reactions are exothermic; the heat 
released during recombination causes higher temperature gradients at the wall resulting in higher 
heat fluxes. 

The skin friction distribution on the body is shown if Fig. 19. It may be noted from this 
figure that the skin friction is zero at 0 = -16.7 deg. indicating that the stagnation point is at this 
location on the body. The distribution also displays two peaks in skin friction values - one at 0 = 
-15.5 deg. and the other at 0 = -17.7 deg. These locations are close to the peak heat flux locations. 
The mass flow in the supersonic jet, after passing through the terminating normal shock, is 
divided into two parts - one going up and the other going down along the body. These accelerat- 
ing flows exiting the stagnation region (between the jet terminating shock and the body) seem to 
cause locally high skin friction. The high skin friction values perhaps lead to the two peaks in the 
heat flux distribution. 

Carlson and Wilmoth predict a peak pressure of 2,900 kPa (approx.) and a peak heat flux 
of 530 MW/m 2 for the shock interference case with the noncatalytic wall (see Fig. 9, Ref. 17). 
The corresponding values from present computations are 2,440 kPa and 43 1 MW/m 2 , and are 
lower than the DSMC computations. The location of the peak pressure on the body also differs - 
the present computations show the peak at -16.7 deg. whereas the DSMC results show the peak at 
around -19 deg. The temperature and Mach number contours in the two flowfields, however, 
appear similar. (Compare Figs. 8a and 8b from Ref. 17 with the present Figs. 15c and 15b). One 
of the significant differences in the two studies is that flow properties in the stream below the inci- 
dent shock are somewhat different in the two cases. In the present computations those conditions 
were obtained by assuming an inviscid flow past a 6 deg. wedge, whereas Carlson and Wilmoth 
computed those conditions using a nonequilibrium DSMC code with no boundary layer on the 
wedge. This would result in differences in the flow properties in that stream which in turn could 
lead to differences in surface pressure and heat fluxes. The present results, however, show con- 
siderable details of the flowfield, particularly near the jet terminating shock and the stagnation 
region. This has been possible due to the very fine grids employed in the present computations. 
Such details are not seen in the DSMC results. 
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Pressure and heat flux distributions were computed for Type IV shock interference in a 
perfect gas (See Prabhu, Ref. 19). Those results were obtained using a highly refined mesh of 
about 235,000 elements, and yielded a peak pressure 2,301 kPa and a peak heat flux of 529 MW/ 
m 2 . Although the impinging shock location in that case was the same as in the present case, the 
normal shock that terminates the jet was not parallel to the body. This is primarily because the 
stand-off distance of the bow shocks in front of the body is much larger for perfect gas than for 
corresponding nonequilibrium flows. Larger stand-off distance in the perfect gas case caused 
wider supersonic jet. The pressure as well as the heat flux distributions for the perfect gas case 
are skewed, suggesting that the impinging shock location did not correspond to the worst stagna- 
tion point heating condition. 

Variations of pressure, density, translational and vibrational temperatures, and mole frac- 
tion of the constituent gases as the flow passes through the jet terminating normal shock and 
reaches the stagnation point (along the Cut A- A shown in Fig. 16) are shown in Figs. 20a - 20d. 
The pressure rises quite sharply as the flow passes through the normal shock after which as the 
flow reaches the surface, the pressure gradually rises to its peak value (see Fig. 20a). The density 
also rises through the normal shock, but rises very sharply through the boundary layer. This is 
because the pressure is essentially constant through the boundary layer whereas the temperature 
drops very rapidly. The translational and vibrational temperature distributions (Fig. 20b) are sim- 
ilar to what was noticed for the undisturbed flow (compare this figure with Fig. 9). Similarly, the 
mole fraction distributions for the fully catalytic wall as well as that for the noncatalytic wall, 
shown in Figs. 20c and 20d, are similar to what was observed for undisturbed flow (see Figs. 8a & 
8b). These similarities are to be expected since in the both the cases the flow is passing through a 
normal shock and reaching a stagnation point at the wall. 

CONCLUDING REMARKS 

The governing equations for chemical and thermal nonequilibrium flow, and some details 
of the solution algorithm which employs a point-implicit, cell-centered, upwind scheme are pre- 
sented. The solution employs Roe’s flux difference splitting method for the inviscid fluxes. The 
method for implementing the scheme on an unstructured meshes is included. Air is treated as a 
reacting gas mixture of 5 chemical species namely 02, N2, 0, NO, and N. Ions and electrons are 
neglected. Two temperatures namely translational-rotational and vibrational are assumed to 
describe the thermodynamic energy. Since the rotational characteristic temperatures for the con- 
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stituent species are small, the translational and rotational energy modes are assumed to be in equi- 
librium, whereas the vibrational energy mode is assumed to be in nonequilibrium. A harmonic 
oscillator model is assumed to characterize the vibrational energy of molecular species. 

The computational domain was discretized using unstructured triangular elements and 
several layers of structured quadrilateral elements on the body. Structured layers of quadrilateral 
elements placed next to the body surface enable capturing large gradients in the flow quantities in 
the boundary layer. The part of mesh with triangular elements is obtained adaptively, and further 
refined for the shock interference case to enable capturing complex flow features. Several layers 
of quadrilateral elements were placed next to the wall to enable capturing the boundary layer fea- 
tures well. The thickness of the first quadrilateral element layer was 2.54xl0~ 6 mm, with a mini- 
mum spacing of 0.056 deg. along the circumference of the body. The final mesh employed had 
155,960 elements - the smallest element had a side length of 0.00254 mm. This may be compared 
with the body radius of 2.54 mm. It has been possible to obtain detailed flow features and heat 
flux and pressure distributions primarily because of this adapted and highly refined meshes. 

Details of the complex flow features associated with Type IV shock interference could be 
seen in the contours plots of the flow variables. The triple points, the shear layers, the supersonic 
jet, and the terminating normal shock may be seen clearly in those plots. For the conditions in the 
example, the results indicate that oxygen is almost fully dissociates behind the bow shock in front 
of the body. For the catalytic wall case the oxygen atoms recombine at the wall. A part of nitro- 
gen also undergoes dissociation; however, the nitrogen atoms recombine at the wall conditions. 
The pressure distribution plot exhibits a small region of high pressure on the body (where the 
supersonic jet impinges). The peak pressure is approximately 13.1 times the pitot pressure in per- 
fect gas. The heat flux distribution exhibits two peaks of almost equal magnitude within a small 
region (also corresponding to the jet impingement region). The peak heat flux is 544 MW/m for 
the catalytic wall and 43 1 MW/m 2 for the noncatalytic wall. The effect of a catalytic wall was 
primarily to increase the heat flux on the body. The corresponding numbers for undisturbed flow 
(with no shock interference) are 20.4 MW/m 2 for the catalytic wall and 15.9 MW/m 2 for the non- 
catalytic wall. 

The adaptive remeshing procedure works satisfactorily for strong shocks, but is unable to 
capture weak shocks and shear layers. It was therefore necessary to refine the regions of shear 
layers and weak shocks manually to enable resolving these flow features. Refining the mesh in 
these regions was necessary for adequately resolving these flow features and obtaining accurate 
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results. The present solution algorithm is only first order accurate, and a higher order solution 
algorithm would perhaps require a less refined mesh. 

Present results for the undisturbed flow case compared well with DSMC computations. 
However, the results for the shock interference case (location and magnitudes of the peak pressure 
and the peak heat flux) differed from the DSMC computations. It is possible that factors includ- 
ing the differences in the post incident shock flow properties and differences in computational 
grids could have contributed to the differences in the results. 
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Table 1: Reaction Rate Constants 


j 

Reaction 

4 

N J 

D j 

1 

02+02 <=> 0+0+02 

23.243e+13 

-1.0 

59,500 

2 

02+N2 <=> 0+0+N2 

7.200e+12 

-1.0 

59,500 

3 

02+0 <=> 0+0+0 

9.000e+13 

-1.0 

59,500 

4 

02+N0 <=> O+O+NO 

3.600e+12 

-1.0 

59,500 

5 

02+N <=> O+O+N 

3.600e+12 

-1.0 

59,500 

6 

N2+02 <=> N+N+02 

1.900e+ll 

-0.5 

113,500 

7 

N2+N2 <=> N+N+N2 

24.700e+ll 

-0.5 

113,500 

8 

N2+0 <=> N+N+O 

1.900e+ll 

-0.5 

113,500 

9 

N2+N0 <=> N+N+NO 

1.900e+ll 

-0.5 

113,500 

10 

N2+N <=> N+N+N 

4.085e+16 

-0.5 

113,500 

11 

NO+02 <=> N+O+02 

3.900e+14 

-1.5 

75,500 

12 

NO+N2 <=> N+0+N2 

23.900e+14 

-1.5 

75,500 

13 

NO+O <=> N+0+0 

27.800e+15 

-1.5 

75,500 

14 

NO+NO <=> N+O+NO 

27.800e+15 

-1.5 

75,500 

15 

NO+N <=> N+O+N 

27.800e+15 

-1.5 

75,500 

16 

NO+O <=> N+02 

23.200e+03 

1.0 

19,700 

17 

N2+0 <=> NO+O 

0.000e+07 

0.0 

38,000 
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Figure 1 . Type IV supersonic jet interference pattern 


y 



Shock r=2.54 mm 

Figure 2. Schematic of the problem 
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(a) Contours of (u/uj (b) Mach number contours 



(c) Translational temperature contours (d) Vibrational temperature contours 


Figure 3. Contours of u-velocity, Mach number, and translational and vibrational temperatures 
for Mach 15 undisturbed flow (fully catalytic wall) 
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(a) 02 mole fraction contours 


(b) O mole fraction contours 



(c) N2 mole fraction contours 



(d) N mole fraction contours 


Figure 4. 02, O, N2, and N Mole fraction contours for Mach 15 undisturbed flow 
(fully catalytic wall) 
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(a) Contours of (u/uj (b) Mach number contours 



(c) Translational temperature contours (d) Vibrational temperature contours 


Figure 5. Contours of u- velocity, Mach number, and translational and vibrational temperatures 
for Mach 15 undisturbed flow (noncatalytic wall) 
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(a) 02 mole fraction contours 


(b) O mole fraction contours 



(c) N2 mole fraction contours 



(d) N mole fraction contours 


Figure 6. 02, O, N2, and N Mole fraction contours for Mach 15 undisturbed flow 
(noncatalytic wall) 
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- 1.30 - 1.25 - 1.20 - 1.15 - 1.10 - 1.05 - 1.00 

x/r 

Figure 8. u-velocity distribution on the centerline for Mach 15 

undisturbed flow (fully catalytic and noncatalytic walls) 
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x/r 

Figure 9. Temperature distributions on the centerline for Mach 15 
undisturbed flow (catalyitc and noncatalytic walls) 
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Figure 10. Pressure distribution on a 0. 1 inch radius cylindrical body 
in Mach 15 flow (catalytic and noncatalytic walls) 



Figure 1 1 . Heat flux distribution on a 0. 1 inch radius cylindrical body 
in Mach 15 flow 
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(b) O mole fraction 


Figure 12. 02 and O mole fraction contours for Mach 15 
shock interference flow; catalytic wall 
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(b) Mach number 

Figure 14. Contours of u-velociy and Mach number for 

Mach 15 shock interference flow; catalytic wall 
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(a) Translational temperature 
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(b) Vibrational temperature 

Figure 15. Translational and vibrational temperature contours 
for Mach 15 shock interference flow; catalytic wall 
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Figure 16. Pressure contours for Mach 15 shock interference flow; 
fully catalytic wall 
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0. deg. 

Figure 17a. Pressure distribution on the body for Mach 15 shock interference flow 
(fully catalytic and noncatalytic walls) 



6, deg. (expanded scale) 

Figure 17b. Pressure distribution near the stagnation point on the body for Mach 
15 shock interference flow (fully catalytic and noncatalytic walls) 
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Figure 18 a. Heat flux distribution on the body for Mach 15 shock interference flow 



Figure 18b. Heat flux distribution near the stagnation point on the body for Mach 15 
shock interference flow 
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Figure 20a. Pressure and desity distributions along Cut A- A for Mach 15 
shock interference flow; noncatalytic wall. 
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Figure 20b. Translational and vibrational tempertaure distributions along Cut A- A 
for Mach 15 shock interference flow; noncatalytic wall 
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d/r 


Figure 20c. Mole fraction distributions along Cut A-A for Mach 15 
shock interference flow; catalytic wall 



Figure 20d. Mole fraction distributions along Cut A-A for Mach 15 
shock interference flow; noncatalytic wall. 
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